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ABSTRACT* 


Opening,  shearing  and  combined  mode  fracture  tests  were  conducted 
with  long  rectangular  strips  of  plexiglas  clamped  on  the  long  edges  and 
containing  centered  and  off-centered  cracks.  The  critical  stress  inten- 
sity factors,  crack  initiation  angles,  and  crack  paths  were  evaluated. 
Fractured  surfaces  were  then  examined  as  to  crack  behavior.  The  maximum 
energy  release  rate  criterion  was  extended  to  problems  with  a large 
degree  of  shearing  mode  present.  This  criterion  was  then  used  to  predict 
successfully  the  initial  and  subsequent  crack  propagation  behavior  pre- 
sented in  the  experimental  work.  Crack  arrest  was  examined  for  this 
rigid  grip  configuration.  Using  in  part  information  obtained  from  these 
studies,  a flexible  fiber  model  was  investigated  to  determine  the  effect 
of  fiber  bending  on  crack  behavior.  Certain  aspects  of  crack  growth  and 
arrest  in  these  idealized  fiber  models  were  explored. 


This  report  is  based  on  a dissertation  submitted  to  the  Graduate  College 
of  Texas  A&M  University  in  partial  fulfillment  of  the  requirements  for 
the  degree  of  Doctor  of  Philosophy  in  Interdisciplinary  Engineering.  The 
author  is  now  at  Haliburton  Services,  Duncan,  Oklahoma. 
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INTRODUCTION 


In  recent  years,  there  has  been  an  increased  interest  in  the  use 
of  composite  materials,  though  the  idea  of  a composite  or  multi- 
phase material  is  centuries  old  [1],  There  are  two  general  groups 
of  composite  materials,  particulate  and  fibrous,  depending  on  the 
form  of  the  added  material  [2],  This  discussion  will  be  limited  to 
fibrous  composites,  consisting  of  stiff,  approximately  parallel  fibers 
in  a softer  matrix  material,  as  in  Figure  1.  The  result  is  an  efficient 
structural  material.  Current  research  is  leading  to  the  use  of 
unusually  strong,  high  modulus  fibers  of  various  materials  with 
polymeric,  ceramic,  or  metallic  matrices,  with  resultant  high  strength- 
to-weight  characteristics  [3,4],  The  most  commonly  used  structural  com- 
posite, however,  consists  of  graphite  or  glass  fibers  in  an  epoxy  matrix. 

These  constituent  materials  are  usually  considered  to  exhibit 
linear  elastic  behavior;  however,  the  matrix  and  composite  may 
exhibit  a significant  amount  of  time  and  temperature  dependent 
mechanical  behavior  in  many  applications.  To  insure  safe  structural 
use  of  these  materials,  the  elastic  behavior  must  be  examined,  then 
extensions  to  viscoelasticity  made.  Much  research  on  polymeric 
composites  has  been  done  in  this  manner  [4,5]. 

Researchers  have  also  found  that  polymeric  composites  are 
nonlinearly  viscoelastic.  Tsai  et  al.  [4],  Lou  and  Schapery  [6], 

Ashton  [7],  and  Schapery  et  al.  [8]  have  shown  these  composites  are  not 
linear  except  at  very  small  stress  levels,  which  may  be  below  struc- 
tural design  limits.  Studies  with  fibrous  composites  indicate  that 
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a sizable  portion  of  this  nonlinearity  is  due  to  internal  crack  growth 
within  the  matrix  or  fiber/matrix  interface  [6,8,9].  This  nonlinearity 
due  to  crack  growth  is  significant,  particularly  during  the  first  one 
or  two  cycles  of  loading  and  unloading  of  the  structure,  where  a 
significant  amount  of  irreversible  damage  occurs.  However,  except  at 
high  stresses  or  temperature,  after  these  initial  cycles  the  response 
does  not  change,  probably  indicating  some  sort  of  crack  arrest  [6,10,11]. 

The  majority  of  the  cracks  originate  as  microscopic  flaws  in  the 
matrix  material  resulting  from  the  manufacturing  of  the  composite 
structure  [12-15].  Upon  loading  these  defects  cause  stress  concentra- 
tions, inducing  crack  initiation  and  growth.  These  cracks  then  pro- 
pagate under  load  (even  below  catastrophic  levels)  in  the  matrix 
through  various  mechanisms  [16-18]. 

The  nature  of  crack  propagation  in  the  matrix  is  very  complex. 
Geometrically,  the  microcrack  in  the  matrix  is  at  some  angle  to  the 
surrounding  fibers  which  are  also  at  some  other  angle  to  the  applied 
loading.  The  problem  of  a crack  in  a low  modulus  material  surrounded 
by  a much  stiffer  material  (microcrack  in  the  matrix  of  an  advanced 
composite)  has  been  investigated  theoretically  [18-21].  This  research 
has  shown  that  either  the  arbitrary  crack  propagates  toward  the 
stiffer  material,  causing  the  stress  intensity  factor  to  drop,  forcing 
the  crack  to  slow  or  arrest,  or  the  crack  may  actually  propagate 
toward  the  center  of  the  matrix  and  then  parallel  to  the  fibers. 

1 he  direction  of  the  fibers  relative  to  an  annlied  tensile  load 
may  have  a significant  affect  on  the  tendency  for  microcrack  arrest 
whether  or  not  the  crack  is  parallel  to  the  fibers:  indeed,  as  we  will 


4 


show  a crack  which  is  initially  aligned  with  the  fibers  could  turn  into 
the  fibers  with  subsequent  arrest  if  the  shear  stress  along  the  fiber 
direction  falls  within  a certain  range. 

Researchers  have  presented  another  possible  mechanism  of  crack 
arrest,  prevalent  even  with  cracks  parallel  to  the  fibers.  Fiber 
bending  may  cause  cracks  to  then  "skip"  across  the  fiber  and  continue 
propagation  in  the  same  plane  as  the  original  crack  [8,22,23]  forming 
a crack  arrest-growth  mechanism. 

The  engineering  theories  of  fracture  mechanics  can  be  used  to 
investigate  the  nature  of  the  crack  arrest  mechanisms  on  a micro- 
structural  scale  in  composite  materials.  The  work  presented  here 
starts  with  a review  of  the  classical  theory  of  fracture  for  a crack 
at  an  arbitrary  orientation  relative  to  the  applied  loads.  Experi- 
mental data  are  then  presented  for  a geometry  similar  to  a cracked 
matrix  surrounded  by  perfectly  rigid  fibers.  Fracture  mechanics, 
experimental  data,  and  fractography  data  are  compared.  A cracked- 
composite  model  is  then  examined  through  a fracture  mechanics  approach 
to  investigate  some  of  the  effects  of  fiber  deformation  in  fibrous 


composite  materials. 


REVIEW  OF  SOME  ASPECTS  OF  FRACTURE  MECHANICS 


General  Development 

The  study  of  cracks,  whether  large  or  small,  and  their  affect  on 
structural  or  material  behavior  is  generally  entitled  fracture  mechanics. 
Fracture  mechanics  is  often  studied  assuming  linear  elastic  behavior, 
neglecting,  for  example,  the  nonlinear  effects  occurring  in  a small  region 
about  the  crack  tip.  However,  there  are  many  investigations  which  treat 
plasticity  and  viscoelastic  effects  [e.g.  24,  25]. 

The  formal  study  of  fracture  mechanics  started  with  Griffith  [26]. 

He  used  a linear,  elastic  stress  solution  provided  by  Inglis  [27] 
for  a flat  plate  under  uniform  tension  with  an  elliptical  hole,  which 
was  then  degenerated  into  a crack.  Without  using  the  near  tip  stress 
field  directly,  Griffith  devised  an  energy-rate  analysis  for  the  pre- 
diction of  initiation  of  crack  growth  in  brittle  materials.  Not  until 
Sneddon  were  the  stress-field  expansions  for  the  crack  tips  analyzed 
[28].  The  general  applicability  of  these  field  equations  were  recog- 
nized later  by  Irwin  [29,30]  and  Williams  [31],  and  expanded  to  a 
general  form  for  the  isotropic  elastic  body. 

Using  the  methods  of  Westergaard  [32],  Irwin  [33]  derived  three 
sets  of  equations  describing  three  basic  types  of  displacement  fields 
as  shown  in  Figure  2.  The  opening  mode  is  associated  with  a local 
displacement  in  which  the  crack  surfaces  move  directly  apart,  normal 
to  the  x-z  plane.  The  anti-symmetric  or  sliding  mode  is  characterized 
by  displacements  in  which  the  crack  surfaces  slide  over  one  another 
perpendicular  to  the  leading  edge  of  the  crack.  The  third  mode,  skew 
symmetric  or  tearing,  is  described  by  the  crack  surfaces  moving  past  one 
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another  parallel  to  the  leading  edge  of  the  crack.  Irwin's  equations 
for  the  local  stresses  and  displacements  (with  r,6  the  local  polar 
coordinates  in  the  x-y  plane)  are: 


Mode  I 


kt  e ,,  , e , 36, 

o = r-pr  cos  77  (1  - sin  77  sin  -r— ) + 

XX  (2irr)1/2  2 2 2 

a = 4-ry  cos  | (1  + sin  | sin  + 

yy  (2Trr)1/2  2 2 2 

KI  , 0 0 30 

°xv  = 1/2  Sln  2 cos  2 COS  T • • • 

y (2nr)  /Z  ^ ^ 

a = v(p  + a ) , a =a  =0 

zz  xx  yy  xz  yz 

^Ir'*'^2  6 20 

u = — (.■—/  cos  2v  + sin  —)  + 

Kt  r m y a 

v = — (— ) sin  (2  - 2v  - cos  — ) + 

w = 0 


(1) 
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Each  of  these  equations  contains  only  the  first  term  of  an  expansion 
about  the  crack  tip  and  is  exact  as  r approaches  zero.  Since  it  has 
been  generally  believed  and  shown  by  others  [34]  that  the  region  about 
the  crack  tip  is  in  a state  of  plane  strain,  the  first  two  sets  of 
equations  are  written  for  this  case. 

The  K' s in  the  above  equations  are  parameters  called  stress 
intensity  factors.  As  these  factors  are  not  dependent  on  the  local 
coordinates  (r,0),  they  control  the  intensity  of  the  stress  field, 
not  the  shape  of  it.  These  factors  depend  on  the  magnitude  of  the 
applied  forces  and  the  geometry  of  the  body  containing  the  crack,  or 
cracks,  including  crack  size  [35,36].  Their  values  can  be  calculated 
from  analyses  that  are  either  three-dimensional,  plane  stress  or  plane 
strain,  depending  on  the  problem  being  considered.  This  leads  to  confu- 
sion since  the  crack  tip  region  is  considered  to  be  in  plane  strain 
[34,  37].  This  paradox  is  due  to  the  attempt  to  investigate  a 
locally  three-dimensional  problem  using  two-dimensional  concepts.  Since 
these  parameters  are  functions  of  the  applied  loads,  their  values  at 
failure  are  considered  as  "critical"  stress  intensity  factors.  These 
values  then  can  be  determined  by  experiments.  Since  the  inception  of 
the  stress-intensity  factor  approach,  many  geometries  have  been 
analyzed  [e.g.,  see  37,38].  Due  to  limitations  of  analytical  tech- 
niques, only  relatively  simple  geometries  have  been  analyzed. 

Numerical  Analysis  in  Fracture  Mechanics 
Cracked  structures  rarely  have  simple  crack  geometries.  This 
has  lead  to  an  increased  amount  of  research  in  the  area  of  numerical 


techniques  needed  to  analyze  cracked  bodies.  Boundary  collocation 
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schemes  and  finite  element  methods  have  been  the  major  approaches. 

The  boundary  collocation  technique  will  be  discussed  first. 

This  technique  starts  with  a stress  function  which  is  used  to  describe 
the  stresses  throughout  a cracked  body  and  satifies  the  boundary 
conditions  at  discrete  points.  This  procedure  offers  some  advantages 
over  other  methods  such  as  little  data  preparation,  great  accuracy, 
and  a smaller  number  of  equations  to  be  solved  [39,40].  The  colloca- 
tion technique  does  not  have  the  flexibility  of  other  methods  in  that 
of  allowing  curved  cracks,  complex  external  boundaries  and  complicated 
loading  conditions. 

Due  to  its  flexibility,  the  finite  element  method  has  been  used 
extensively  in  fracture  analysis.  The  finite  element  method  itself 
is  well  documented  and  will  not  be  discussed  here  [41].  However, 
a brief  review  of  the  various  approaches  to  fracture  mechanics  with 
the  use  of  finite  elements  will  be  presented. 

One  approach  is  called  the  displacement  method.  Displacements 
are  obtained  from  a model  with  a coarse  outer  grid  and  a very  fine 
grid  near  the  crack  tip.  By  using  the  calculated  nodal  displacements 
along  the  crack  face  and  the  displacement  equations  about  a crack  tip, 
described  previously,  a stress  intensity  factor  can  be  estimated  from 
each  nodal  point.  Since  the  displacements  are  not  exact,  a curve  of 
estimated  K's  versus  the  location  of  the  nodal  point  with  respect  to 
the  crack  tip  is  obtained,  rather  than  a single  value  for  K.  Chan 
et  al.  [42]  have  found  that  the  curve  appeared  to  be  linear 
some  distance  from  the  crack  tip,  the  distance  decreasing  with 
increased  grid  refinement.  A tangential  extrapolation  on  the  curve 


to  the  crack  tip  was  used  to  estimate  K.  The  answers  obtained  in 
this  manner  can  be  reasonably  accurate  (<5%  error)  but  only  with 
very  refined  meshes  (elements  on  the  order  of  element  area/crack 
length  squared,  (3.2  x 10  ) [42].  Obviously,  this  method  is  not 
very  precise  in  execution  and  very  costly  by  requiring  such  fine 
meshes. 

Others  [43,44]  have  added  higher  order  terms  to  the  local  crack 
tip  displacement  equations  with  some  success.  Using  the  logic  that 
as  the  calculated  displacements  increase  in  accuracy  and  the  local 
displacement  equations  decrease  in  accuracy  as  one  moves  from  the 
crack  tip,  a suitable  compromise  might  be  achieved.  Plotting  as 
before,  a maximum  K can  be  found  and  can  also  be  fairly  accurate 
(<5%  error).  Mesh  refinement  is  still  a problem  as  is  the  preciseness 
of  execution. 

A novel  approach,  also  using  displacements,  has  recently  been 
proposed  [45].  This  method  is  based  on  using  the  calculated  nodal 
displacements  to  plot  the  deformed  shape  of  the  crack  and  compare  it 
to  the  ellipse  predicted  by  an  elastic  fracture  mechanics  analysis  for 
an  elliptical  crack.  This  method  is  also  limited  in  requiring  some 
mesh  refinement  and  is  only  applicable  to  opening  mode  cracks. 

An  approach  similar  to  the  displacement  method  is  based  on  using 
calculated  stresses.  The  stresses  can  be  handled  as  were  the  displace- 
ments, by  comparing  with  the  stress  singularity  equations  [42].  Since 
stresses  are  less  accurate  than  displacements  for  codes  using  the 
method  of  direct  stiffness,  this  approach  suffers  an  additional  source 
of  error  [46],  This  method  and  the  displacement  method,  however,  are 
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theoretically  applicable  to  three  dimensional  problems. 

Another  approach  using  conventional  finite  element  codes  employs 
line  integrals  such  as  Rices' s .7  integral  [471.  This  line  integral  is 

3u . 

J = / (Wdy  - T.  ~ ds)  (4) 

Y i dx 

where  v is  an  arbitrary  contour  surrounding  the  crack  tip,  W is  the 
strain  energy  density,  T are  the  surface  tractions,  the  ^ are  the 
displacements,  and  ds  is  an  element  of  the  arc  length  along  . This 
expression  can  also  be  related  to  the  stress  intensity  factors  for 
plane  strain  problems  by 

J = (1  - v2)(Kt2  + K:i2)/E  (5) 

where  E is  Young's  modulus.  This  approach  does  not  require  as  fine 
a mesh  as  the  previous  methods,  but  the  contour  must  be  adjusted  to 
minimize  the  error  [46].  The  integral  is  theoretically  independent 
of  contour  path  but  is  path  dependent  in  practice.  The  method  also 
suffers  in  that  only  one  mode  of  fracture  can  be  analyzed  at  one  time, 
and  the  crack  surfaces  can  not  be  loaded  [43,48,49]. 

Another  contour  integral  has  been  developed  which  can  solve 
combined  Mode  I and  II  problems  (where  both  and  are  calculated), 
even  for  anisotropic  solids  [50].  Using  a form  of  Betti's  reciprocal 
work  relation  and  comparing  contours,  a path  independent  integral  can 
be  derived  that  is  fairly  accurate. 

A different  approach  to  fracture  problems  using  finite  elements  is 
based  on  the  strain  energy  release  rate,  G,  where 


which  is  the  change  in  strain  energy  dW  for  a change  in  crack  surface 
area  dA.  This  quantity  can  be  related  to  other  common  terms  as 


G = J = (1-v2)(Kx2  + Kit2)/E 


(7) 


for  plane  strain  and  more  generally  as  [30] 


r ,,  2WV  2 . „ 2 III 

G=(i-v)(KI  + Kn  +-(Tr^y)/E 


(8) 


By  obtaining  two  finite  element  solutions  of  the  same  body  with  two 
slightly  different  crack  lengths,  G can  be  calculated.  This  method 
has  been  sucessfully  developed  and  used  by  many  researchers  [51-54]. 
Good  accuracy  may  be  realized  with  course  meshes.  This  is  due  to  the 
accuracy  of  the  finite  element  method  in  calculating  strain  energy  and 
the  cancelling  of  errors  when  taking  the  difference  in  strain  energy 
for  two  different  lengths  of  colinear  cracks.  In  order  for  the  errors 
to  cancel,  the  local  mesh  must  be  identical  for  both  problems.  As  with 
the  J-integral,  this  approach  permits  the  determination  of  individual 
stress  intensity  factors  when  only  one  mode  exists. 

Another  strain  energy  approach  is  to  relate  i t to  the  rate  of 
change  of  structural  compliance,  C,  or  the  inverse  spring  constant, 
with  crack  extension  [55]: 


r . P‘!  dC 
G dX 


(9) 


where  P is  the  applied  load.  Pook  and  Dixon  [55]  calculated  the  compli- 
ance of  a cracked  body  for  a number  of  crack  lengths  by  the  use  of  finite 


14 

element  codes  and  then  used  (9)  to  calculate  G.  This  approach  is 
comparable  to  the  previous  one,  but  the  strain  energy  is  now  calcu- 
lated as  one-half  the  summation  of  the  product  of  applied  forces  times 
their  respective  change  in  displacements.  The  accuracy  is  good,  but, 
again,  only  one  mode  at  a time  can  be  analyzed  for  K. 

Since  the  form  of  the  elastic  crack  tip  singularity  is  known 
through  equations  (1-3) , special  crack  tip  elements  may  be  constructed 
which  have  the  singularity  condition  built  into  the  local  displace- 
ment pattern.  Use  of  such  an  element  eliminates  a need  for  fine  meshes 
near  the  crack  tip  and  avoids  the  questionable  convergence  of  the 
finite  element  method  at  the  crack  tip.  There  are  many  such  methods 
presented  in  the  literature  and  they  can  be  divided  into  special  ele- 
ments surrounding  the  crack  tip,  special  elements  connecting  at  the 
crack  tip,  or  modified  conventional  elemencs. 

The  elements  of  the  first  type  were  developed  initially  by  Wilson 
[56]  and  bv  Hilton  [57],  where  displacements  generated  from  William's 
displacement  equations  [31]  were  used  as  boundary  conditions  for  con- 
ventional elements.  These  fracture  elements  which  are  circular, 
contain  a centered  crack  tip  extending  radially  to  the  outer  edge. 
Wilson's  element  is  for  Mode  III  problems  and  Hilton's  for  Mode  I. 

The  displacements  in  the  elements  are  exact  as  r-H),  thus  a balance 
must  be  achieved  between  the  abilities  of  the  special  elements  and 
conventional  elements  to  model  the  overall  displacement  field. 

Wilson  extended  his  work  to  plane  problems  [46,58]  including  higher 
order  terms  in  the  displacement  equations,  with  some  improvement  in 
accuracy.  Byskov  [59]  has  a similar  triangular  element  based  on 
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displacements  derived  from  the  complex  stress  function  of  Muskheli shvi li . 
In  all  the  above  elements,  there  is  no  compatibility  between  the 
singularity  elements  and  conventional  elements  along  their  common 
sides.  Wilson's  element  can  have  an  increasing  number  of  conventional 
elements,  leading  to  a pseudo-convergence,  but  Byskov's  element  is 
coupled  with  only  three  constant  strain  triangular  elements  [46]. 

Others  have  also  developed  similar  approaches  using  multiple  special 
regions  [60,61]. 

Tong  et  al.  [62]  have  derived  a crack  tip  element  based  on 
the  hybrid  element  concept  that  is  compatible  with  conventional 
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elements  [63].  Using  assumed  stress  distributions  based  on  the  r 
type  streps  field,  a stiffness  matrix  was  derived  for  an  octagonal 
element  (may  be  rectangular).  Linear  displacement  is  assumed  at  the 
boundary  so  that  the  resulting  displacements  will  then  be  compatible 
with  the  displacements  from  conventional  elements  at  their  common  sides. 
Due  to  the  nature  of  the  derivation,  the  stress  intensity  factors  are 
an  immediate  result.  This  element,  for  Mode  I and  II,  is  readily  used 
with  coarse  meshes  and  is  easily  implemented  with  other  finite 
element  codes  [64].  Others  have  also  implemented  similar  procedures 
to  derive  other  elements  [65,66], 

Researchers  have  also  developed  triangular  elements  connecting 
at  the  crack  tip,  notably  Wilson  [58]  and  Tracey  [67].  These  elements 
have  a series  of  triangles  centered  at  a crack  tip  much  like  a cut 
up  pie.  In  Wilson's  formulation,  the  displacements  vary  linearly 
with  respect  to  the  angle  about  the  crack  tip  and  are  proportional 
to  the  square  root  of  the  radius.  Displacement  is  continuous  at  the 
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Interface  of  the  crack  elements  but  are  continuous  with  the  conventional 
elements  at  only  the  nodes.  This  lack  of  compatibility  occurs  because 
the  displacements  on  the  radial  side  of  the  crack  tip  element  do  not 
vary  in  the  same  way  as  with  conventional  elements.  As  more  elements 
are  used  in  the  overall  singular  crack  tip  element,  the  displacement 
over  the  radial  edge  of  the  element  approaches  a linear  relation  and 
continuity  with  surrounding  elements  is  approached.  Tracey's  element 
allows  the  displacements  to  vary  linearly  at  the  radial  edge  assuring 
continuity  with  the  conventional  elements.  Neither  formulation  is 
easily  implemented. 

The  third  approach  to  modeling  stress  singularities  is  to  modify 

existing  finite  elements  close  to  the  crack  tip.  A recent  method 

utilizes  quadratic  isoparametric  elements.  The  modified  element's 

midside  nodes  on  the  sides  near  the  crack  tip  are  moved  to  the  quarter 

points  near  the  crack  tip  [68-70] . The  moving  of  the  node  causes  a 
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singularity  in  the  form  of  r at  the  crack  tip  node  for  strains, 
thus  a special  element  is  formed.  Either  the  strain  energy  release 
method  or  the  displacement  technique  can  then  be  used  to  calculate 
the  stress  intensity  factors.  This  method  seems  to  be  easily  used 
and  fairly  accurate,  but  is  not  quite  as  good  as  hybrid  elements  [67]. 
Three  dimensional  possibilities  are  also  being  investigated  [71]. 

There  are  also  a large  number  of  other  modified  finite  element 
techniques.  Some  are  based  on  Buekner's  weighting  functions  [72,73], 
others  use  transition  functions  [74],  substructuring  techniques 
[75,76],  or  elements  with  artificially  low  moduli  [77].  Plasticity 
effects  have  been  considered  [57,78]  and  anisotropic  fracture  is 


being  investigated  by  finite  element  techniques  [79-81]. 

The  previous  discussion  covered  the  analysis  of  stationary  cracks. 
Also  of  great  interest  is  the  criterion  for  initiation  of  crack  growth 
and  the  direction  of  propagation. 

Mixed  Mode  Fracture  Mechanics 

According  to  Griffith,  "the  crack  will  grow  in  the  direction  along 
which  the  elastic  energy  release  per  unit  crack  extension  will  be 
maximum  and  the  crack  will  start  to  grow  when  this  energy  reaches  to 
a critical  value"  [82],  At  first,  analytical  complexities  prevented 
the  analysis  of  non-colinear  crack  extension  problems,  so  only  simple 
Mode  I problems  were  investigated.  Since  real  structural  problems 
often  involve  combinations  of  two  or  three  modes,  one  must  be  able  to 
account  for  the  situation  in  which  the  crack  growth  direction  is  not 
in  the  original  plane  of  the  crack.  A summary  of  the  various  mixed 
mode  criteria  will  be  presented  here. 

One  of  the  first  studies  in  this  area  was  conducted  by  Erdogan  and 
Sih  [82]  in  the  investigation  of  non-colinear  crack  growth  in  two- 
dimensional  geometries.  The  investigation  supported  a criterion  stating 
that : 

1)  Crack  extension  starts  at  its  tip  in  a radial  direction 

2)  Crack  extension  starts  in  a plane  perpendicular  to  the  direc- 
tion of  greatest  local  tension. 

Using  the  singular  portion  of  the  equations  described  earlier  for  the 
stresses  at  a crack  tip,  but  in  polar  coordinates. 
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the  angle  of  propagation  can  be  found  by  maximizing  o00  with  respect 
to  0 as  in 


,1/2  9 CT00 
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(11) 


where  r is  factored  out  of  the  equation.  This  results  in 


cos  — [ sin0  + Kjj (3cos0  -1)]  = 0 


(12) 


which  gives 


0 = + 71 


and 


Kj  sin  0 + ^(3  cos  0 -1)  = 0 (13) 

as  solutions.  The  first  solution  corresponds  to  the  free  surface  condi- 
tions of  the  crack  and  is  discarded;  (13)  yields  the  maximum  o Know- 

00  . 

ing  the  particular  stress  intensity  factors  of  the  problems,  prior  to 
growth,  the  crack  trajectory  can  be  determined.  This  criterion  was 
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compared  to  experiments  conducted  with  plexiglas  sheets  with  various 
crack  geometries  and  loading  conditions.  Most  widely  used  was  a large 
sheet  containing  a centered  crack  inclined  to  the  uniform  tension  at 
the  edge  of  the  sheet.  The  results  compared  fairly  well  with  the 
trajectory  criteria,  allowing  for  the  large  amount  of  data  scatter. 

An  important  comment  was  made  [81]  in  that  the  greatest  energy  release 
occurs  when  the  crack  propagates  approximately  perpendicular  to  the 
direction  of  maximum  tension,  indicating  another  possible  criterion. 

The  same  problem  was  reinvestigated  by  J.G.  Williams  and  P.D. 

Ewing  [83]  with  a slightly  improved  criterion.  Inclusion  of  the  next 
higher  order  term  in  the  equation  for  o_.  improved  correlation  with 

U D 

the  previous  experimental  evidence  and  also  with  new  data  gathered  by 
the  investigators.  Again  a large  amount  of  data  scatter  was  present 
and  the  same  inclined  crack  geometry  with  plexiglas  was  tested.  A 
correction  was  made  by  Finnie  and  Saith  [84]  to  the  equations  Williams 
and  Ewing  used  which  also  improved  correlation  with  data.  The  added 
term  approach  has  forced  the  necessity  of  evaluating  the  expressions 
at  a certain  critical  value  of  the  radial  distance  from  the  tip,  r, 
while  previously  r had  been  factored  out  as  in  equation  (11).  The 
value  of  r is  then  adjusted  to  fit  experimental  data.  In  further 
comments,  Williams  and  Ewing  pointed  out  that  even  higher  order  terras 
may  be  needed  for  bending  experiments  [85].  This  was  followed  up  with 
tension  and  bending  tests  with  inclined  edge  cracks  and  comparison  with 
the  higher  order  term  approach  met  with  some  success  [86].  For  com- 
pleteness, Sih  and  Kipp  [87]  conducted  an  analysis  examining  the  stresses 
found  from  the  exact  solution  computed  from  a degenerate  elliptical 
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cavity.  The  use  of  these  equations  is  also  dependent  on  evaluation 
with  a certain  value  of  r.  This  criterion  compared  with  experimental 
data  with  about  the  same  success  as  the  previous  studies  [87].  Some 
of  these  maximum  stress  criteria  are  presented  in  Figure  3 for  only 
one  value  of  r. 

A minimum  strain  energy  criterion  has  been  proposed  by  Sih 
[87-90].  The  work  is  an  extension  of  some  ideas  presented  earlier 
[82].  The  change  in  strain  energy  density  with  polar  angle  at  the 
crack  tip  was  derived  as  a homogeneous  quadratic  form  of  and 
This  new  criterion  is  specified  as  follows: 

1)  The  crack  will  propagate  in  the  direction  of  minimum  strain 
energy  density  at  the  crack  tip 

2)  The  onset  of  growth  will  be  controlled  by  a critical  value 
of  the  minimum  strain  energy  density. 

Using  only  the  singular  portion  of  the  stress  equations  (1-3), 
the  strain  energy  density  of  a infinitesimal  cube  is  calculated, 
resulting  in: 


§ ' ^2  <“ll  KI2  + 2 a12  Vll  + “22  hi'  + “33  KIII2>  <14> 

ri: 

where  the  coefficients  a^  are  functions  of  the  polar  angle  about  the 
crack  tip.  Young's  modulus  and  Poisson's  ratio.  This  three  dimensional 
criterion  is  written  as: 


5 ‘ I?  <“ll  KI2  + 2 “12  Vll  + “22  KII2  + “33  KIII2>  (15) 


where  S is  proportional  to  the  strain  energy  density.  This  equation 
is  examined  to  find  its  minimum  with  respect  to  the  polar  angle  and  the 


Figure  3.  Crack  propagation  angle  vs  stress  intensity  factor  ratio  for  various 
maximum  stress  criteria. 
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value  of  the  angle  producing  the  minimum  is  the  direction  of  propagation. 

This  criterion  has  been  compared  with  experimental  evidence  in 
tension  and  compression  [86,90-92]  and  has  been  found  to  have  little 
improvement  over  the  previous  maximum  stress  approach.  The  S criterion 
has  been  questioned  for  its  dependence  on  Poisson's  ratio,  v;  it  has 
failed  to  predict  accurately  crack  growth  direction  in  cross-rolled 
beryllium  where  the  value  of  Poisson's  ratio  approached  zero  [91]. 
in  this  case,  the  maximum  stress  criterion  was  significantly  better. 

The  minimum  strain  energy  density  criterion  has  also  shown  poor  correla- 
tion with  cone  fractures  in  fused  silica  where  the  material  is  nearly 
incompressible  (v=.5)  [92].  It  also  has  not  been  proven  that  the  criti- 
cal value  of  S is  a material  constant  as  it  should  be  [90], 

In  order  to  improve  the  correlation  with  data,  a more  general 
and  physically  realistic  model  has  been  proposed  by  Sih  [87,89,93]. 

The  crack  is  approximated  by  a narrow  ellipse.  The  exact  energy 
density  based  on  linear  elastic  theory  is  calculated  with  the 
complete  series  expansion  of  the  stresses  about  the  ellipse  crack  tip 
[87].  This  modified  criterion  is  then  dependent  upon  the  radius  of 
curvature  of  the  ellipse,  p,  and  the  distance  from  the  crack  tip  at  which 
the  modified  S is  evaluated.  Sih  explains  that  this  critical  value  of  r 
depends  on  the  crack  geometry  and  loading  [89].  By  adjusting  P and  r, 

Sih  has  demonstrated  improved  correlation  with  data  [87,92]. 

Francis  and  Ko  [94]  have  examined  the  minimum  strain  energy  den- 
sity criterion  in  terms  of  the  stresses  about  an  elliptical  crack. 

They  derived  three  forms  based  on  the  use  of  plane  strain,  plane  stress, 
and  deviatoric  plane  stress  equations.  As  shown  in  Figure  4,  for  one 


Figure  4.  Crack  propagation  angle  vs  stress  intensity  factor  ratio  for  various 
minimum  strain  energy  density  criteria  evaluated  for  v=.3. 
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value  of  Poisson's  ratio,  there  is  some  difference  between  the  various 
forms  of  the  S criterion. 

Another  approach  to  the  combined  mode  problem  is  to  go  back  to 
Griffith's  quote  given  earlier  and  to  analyze  the  energy  release  as  a 
crack  grows  a very  small  amount.  Two  investigations  have  been  made 
using  the  advanced  analytical  tools  developed  since  Griffith  to  examine 
the  combined  Mode  I and  II  problem.  Hussain  et  al . [95]  started  with 
a line  crack  which  has  an  angular  extension  at  one  end.  Considering  a 
virtual  growth  of  this  extension  and  using  the  J-integral  [47]  and 
equating  the  integral  to  strain  energy  release  rate,  a maximum  of  the 

release  was  found  with  respect  to  the  angle  of  the  extension.  The 
authors  use  various  analytical  techniques,  such  as  mapping  functions, 
complex  potentials  and  limit  procedures. 

Due  to  some  questionable  assumptions  and  failure  to  conform  with 
the  known  stress  and  displacement  fields  in  the  limit  of  a straight 
crack  that  were  pointed  out  by  Knauss  and  Palaniswamy  [92],  an  improved 
analysis  was  conducted  in  [92].  Using  a complex  potential  for  the  stress 
analysis  of  a crack,  and  a crack  with  various  angled  increments  at  one 
end,  the  strain  energies  of  various  configurations  were  evaluated. 

Then  the  release  rate  was  calculated  and  a maximum  found  with  respect 
to  the  angle.  The  length  and  angle  of  the  increments  were  varied  to 
find  the  maximum.  This  analysis  was  correlated  with  experimental 
evidence  with  a good  deal  of  agreement  [92];  however,  the  length  of  the 
increment  had  to  be  adjusted  to  fit  the  data.  The  physical  basis  for 
using  other  than  a vanishingly  small  amount  of  crack  growth  is  not 
clear.  The  analysis  was  not  carried  out  for  the  total  range  of  the 
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problem,  leaving  undone  the  Important  case  where  shear  predominates 
in  the  combined  problem;  see  Figure  19.  These  maximum  energy  release 
criteria  are  compared  in  Figure  5. 

Other  researchers  have  used  finite  element  codes  to  examine  the 
combined  mode  problem.  These  have  followed  the  maximum  strain  energy 
release  rate  and  maximum  stress  criteria  for  plane  fracture  problems. 
Coughlan  and  Barr  analyzed  the  infinite  sheet-inclined  crack  problem 
using  the  energy  release  criterion  [96]  and  found  results  that  fall 
between  the  maximum  stress  criterion  and  the  strain  energy  density 
criterion.  As  seen  in  Figure  5 the  results  agree  well  with  the 
maximum  energy  release  rate  criterion  in  some  portions  of  the  plot. 
Unfortunately,  the  finite  element  analysis  was  also  not  followed 
through  in  problems  where  large  amounts  of  shear  are  present.  Marinshaw 
and  Lindsey  [97]  performed  an  investigation  of  the  edge-notched  biaxial 
strip  with  experiments  with  solid  rocket  propellant  and  case  liner 
materials.  Predictions  were  made  of  the  crack  angle  using  a finite 
element  code  to  calculate  the  maximum  normal  stresses  about  the  crack 
tip.  Incorrect  propagation  angles  were  predicted  by  the  classical 
sharp  tip  solutions  of  linear  theory  and  finite  elements,  due  at  least 
in  part  to  the  blunt  nature  of  the  actual  tip.  Use  of  an  artificially 
blunt  crack  tip  in  the  finite  element  analysis  produced  a favorable 
correlation.  It  may  be  conjectured  that  stress  equations  for  ellipses 
could  provide  equally  good  correlation.  Helen  and  Blackburn  [98]  have 
also  conducted  a finite  element  investigation  using  the  energy  release 
rate  method.  Their  results  are  very  poor  for  cases  where  shear  is  more 
significant  than  tension  [98]. 
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Other  analyses  are  being  proposed  for  mixed  mode  fracture  mechanics. 
Silby  and  Cardew  have  [99]  developed  a theory  comparing  the  stress 
intensity  factors  at  both  the  initial  crack  tip  and  the  incremented 
crack  tip  by  means  of  quadratures  [100],  This  three-dimensional 
formulation  compares  well  in  the  two  dimensional  case  with  certain 
experimental  data  [99].  Morozov  and  Fridman  [101]  have  made  a suggestion 
that  cracks  generate  and  propagate  along  a geodesic  line  (shortest 
distance  between  two  points  on  a surface)  passing  through  the  point 
of  initiation.  An  approximate  variational  scheme  was  proposed  based 
on  the  minimization  of  the  difference  between  the  surface  energy  and 
the  strain  energy  of  the  body.  Lindsay  [102]  has  proposed  a hybrid 
criterion  using  the  crack  growth  orientation  as  found  from  the  maximum 
stress  criterion  and  a failure  law  based  on  the  stress  invariants  near 
the  crack  tip. 

On  the  basis  of  the  literature  review,  it  is  the  belief  of  this 
author  that,  for  sufficiently  small  regions  of  nonlinearity  at  the 
crack  tip,  the  maximum  energy  release  rate  theory  governs  crack 
growth  at  all  stages  of  growth.  The  maximum  stress  criterion  represents 
an  approximation  which  improves  as  the  crack  propagates  and  aligns 
itself  perpendicular  to  the  greatest  tension.  An  attempt  was 
made  to  prove  that  the  two  criteria  are  equivalent  [103]. 

For  solving  engineering  problems,  a surface  relating  the  stress 
intensity  factors  (K^,  K^,  Kj^)  can  serve  as  an  initiation  cri- 
terion. Figure  6 shows  such  a criterion  for  two  dimensional  problems 
[82,  88,  94,  95].  The  energy  release  rate  criterion  fits  experi- 
mental data  quite  well  as  shown  in  a later  section  of  this  paper. 
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Figure  6.  Interaction  of  K and  K for  various 
criteria. 
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However,  it  still  remains  to  calculate  the  energy  release  rate  for  the 
complex  three  dimensional  piOblem. 
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FRACTURE  TESTS 
Motivation 

An  experimental  program  was  designed  to  provide  fracture  data  to 
investigate  the  nature  of  crack  growth  between  fibers  in  a composite 
material.  This  data  was  also  to  be  compared  to  predictions  made  by 
the  various  crack  trajectory  theories  to  see  which  criterion  fit  the 
experimental  data  best.  Thus  the  experimental  model  should  be  indica- 
tive of  a cracked  matrix  restrained  by  relatively  rigid  fibers  and 
should  be  a new  test  of  combined  mode  fracture  theories.  While  micro- 
cracks in  composites  are  truly  three  dimensional,  a biaxial  strip 
with  initial  cracks  parallel  to  the  long  sides  was  considered  to  be  a 
suitable  model  to  begin  investigation  of  microcracking  phenomena. 

Material  Selection 

The  microcracking  behavior  described  earlier  is  found  in  most  com- 
mon fiber  reinforced  composites:  graphite  or  glass  fibers  in  a plastic 

matrix.  These  plastics  usually  are  epoxy  systems  that  have  been  especi- 
ally formulated  for  use  in  composites  [1,2]. 

In  the  interests  of  being  able  to  compare  data  with  the  litera- 
ture and  ease  in  use,  plexiglas  was  employed  instead  of  an  epoxy. 

Overall  mechanical  properties  of  the  two  materials  are  quite  similar 
[104,105].  The  fracture  behavior  could  also  be  expected  to  be 
similar  even  in  the  matter  of  crazing  [106-113]  (crazing  is  the  lessen- 
ing in  density  due  to  microvoid  growth,  usually  found  around  a crack 
tip).  After  surveying  the  literature,  it  was  concluded  that 
crazing  or  ductility  might  be  a factor  in  determining  the  overall 
fracture  direction  [114-115].  Therefore,  in  order  to  compare 
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experimental  evidence  on  a new  geometry  with  the  wealth  of  data 
generated  with  plexiglas  [82,83,85,86]  the  same  material  was  chosen. 

The  ease  in  producing  readily  visible  cracks  was  also  an  important 
factor  in  choosing  plexiglas. 

Samples  and  Equipment 

GM  grade  plexiglas  sheets  as  specified  by  the  Rohm  and  Haas 
Company  [117]  were  ordered  in  three  nominal  thicknesses  (0.03,  0.06, 
0.125  inches).  All  the  sheets  needed  were  ordered  at  one  time  to  reduce 
any  batch  variability.  The  material  was  stored  in  an  area  where  the 
temperature  and  humidity  were  the  same  as  the  testing  area  (75°F  at 
50%  relative  humidity). 

The  sheets,  which  arrived  with  protective  masking  tape,  were  pre- 
pared for  a cutting  to  the  required  dimensions,  by  drawing  the  sample's 
dimensions  directly  on  the  tape.  All  sample  dimensions  were  drawn 
slightly  oversize  to  allow  for  final  trimming.  To  compensate  for  any 
induced  anisotropy  caused  by  manufacturing,  half  of  the  samples  of 
one  sheet  were  cut  perpendicular  to  the  other  half.  Using  a tungsten 
carbide  band  saw  blade,  with  a slow  feed  rate  and  band  saw  speed,  the 
samples  were  rough  cut  from  all  the  sheets.  After  removing  the  tape, 
several  coupons  were  mounted  in  a vertical  mill  at  one  time  and  machined 
to  the  final  dimensions  using  a fly  cutter  tool.  Inspection  of  the 
edges  showed  little  damage  if  care  was  taken  to  prevent  the  samples 
from  vibrating.  The  thinnest  sheets  (0.03  inches)  showed  a tendency 
to  edge  crack  upon  cutting;  samples  with  edge  cracks  were  rejected. 

Constant  crosshead  rate  tests  were  conducted  using  a few  special 
samples  taken  randomly  from  the  three  sized  sheets.  These  samples 
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(10  x 1.25  inches)  were  cut  as  the  above.  Micro-measurements'  strain 
gages,  type  EA-06-125BZ-350,  were  bonded  to  the  sample's  center  both 
■ along  the  axis  and  perpendicular  to  the  axis  of  the  sample.  Mirror 

image  gages  were  applied  to  the  opposite  side  of  the  coupon.  Micro- 
measurements' M-bond  AE-15  adhesive  was  used  and  the  gage-adhesive 
was  cured  under  slight  pressure  for  four  hours  at  140°F,  well  below 
any  temperature  that  could  cause  permanent  dimensional,  physical  or 
chemical  changes  to  this  material.  Strain  gage  measurements  were 
recorded  using  a B&F  Instruments'  Strain  Gage  Acquisition  System  and 
a modified  Hewlett-Packard  561  Digital  Printer.  This  system  operated 

I 

at  the  rate  of  two  channels  per  second,  with  capability  up  to  ten 
channels  of  strain  gage  signal  conditioning  with  wire  length  compen- 
sation. 

The  samples  were  placed  with  their  greatest  length  in  the  loading 
direction  and  clamped  in  the  standard  wedge-action  grips  of  the  constant 
crosshead  rate  Instron  Testing  Machine.  Using  a rate  of  0.02  inches/ 
minute  and  at  the  temperature  and  humidity  described  previously  (slight 
variations  in  these  conditions  were  recorded)  the  samples  were  pulled 
apart  until  sudden  failure  occured  at  the  center,  under  the  strain 
gages.  During  the  test  the  applied  crosshead  displacement,  resulting 
i load,  the  two  axial  strains  (both  sides)  and  the  transverse  strains 

(both  sides)  were  recorded. 

In  reducing  this  data,  the  opposing  strains  were  averaged.  The 
transverse  strains  were  adjusted  for  transverse  wire  effects  in  the 
gages  due  to  the  applied  displacement,  by  equations  provided  by  Micro- 
measurements [118].  Using  the  stress  (the  resultant  load/original  area) 
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and  the  strain  (averaged  axial)  from  these  tests,  a stress/strain 
curve  was  plotted.  This  is  shown  in  Figure  7.  The  stress/strain 
curve  is  linear  below  the  level  of  applied  stresses  used  in  the  other 
experimental  work.  The  calculated  Young's  Modulus  (o/e  at  a point) 
was  433,000  psi,  comparable  to  the  manufacturer's  value  of  450,000 
psi . The  Poisson's  ratio  was  found  to  agree  with  the  manufacturer's 
value  of  0.35.  The  calculated  values  were  used  in  the  work  in  this 
report . 

With  these  preliminary  tests  completed,  a review  of  the  litera- 
ture [119-122]  led  to  the  biaxial  strip  size  of  8 x 1 inches;  these 
dimensions  would  cause  an  intial  crack  to  be  in  the  biaxial  stress 
field  (far  enough  from  the  ends  of  the  strip).  The  samples  were  cut 
in  the  same  manner  as  the  previous  discussion.  Based  on  Mueller' s 
work  [121]  crack  lengths  of  one  and  one-half  inches  were  decided  upon; 
this  would  let  the  test  configuration  approach  the  geometry  of  a crack 
in  an  infinite  strip. 

The  cracks  were  cut  parallel  to  the  long  sides  of  the  strips, 
as  in  Figure  8.  An  initial  crack  or  slit  was  machined  with  the  use 
of  a very  fine  circular  saw  blade  mounted  in  a vertical  mill.  The 
resulting  slit  was  nominally  0.005  inches  thick.  This  machining 
produced  a fine  slit  with  a small  curvature  at  each  end,  the  curvature 
being  more  pronounced  with  the  thicker  samples.  (Due  to  the  circular 
nature  of  the  blade  and  mounting  of  it,  the  slit  could  not  have  ends 
perpendicular  to  the  faces  of  the  sheet).  A natural  crack  tip  was 
formed  by  first  applying  in-plane  compression  on  the  long  sides  of 
the  strip  ahead  of  the  crack  (with  a vise)  and  second  by  wedging  open 


the  crack  (with  a razor  blade)  to  produce  controllable  colinear  crack 
growth.  This  extension  was  for  a distance  greater  than  the  sheet 
thickness.  This  crack  tip  had  to  be  perpendicular  to  both  sheet 
faces  as  examined  by  the  naked  eye,  or  the  sample  was  rejected.  The 
process  was  repeated  for  the  other  end  of  the  slit. 

The  test  apparatus  was  designed  so  as  to  test  the  strips  in  ten- 
sion, simple  shear,  or  in  a combination  of  the  two.  As  seen  in  Figure 
9 two  plastic  strips  were  bonded  perpendicular  to  stainless  steel  bars. 
The  steel  bars  were  designed  to  be  extremely  rigid  in  respect  to  the 
plexiglas.  The  outer  two  bars  were  wider  than  the  inner  bar  to 
facilitate  the  shear  experiment.  The  perpendicular  bonding  was  assured 
by  a special  bonding  fixture  as  seen  in  Figure  10.  The  bar  was  placed 
in  the  bottom  of  the  fixture,  and  the  inner  moveable  plates  supported 
the  plastic  perpendicular  to  the  bar.  With  the  use  of  a micrometer, 
the  moveable  plates  were  adjusted  until  the  biaxial  strip  was  centered 
perpendicular  to  the  bar.  The  moveable  plates  were  locked  into  place, 
the  two  sides  of  the  fixture  loosened,  the  biaxial  strip  removed,  adhe- 
sive applied,  the  system  reassembled  and  left  until  bonding  was  complete 

Two  strips  of  identical  dimensions  were  bonded  to  two  bars 
and  then  the  bars  placed  so  as  to  have  the  strips  parallel  to  the 
table's  surface.  The  smaller  middle  bar  was  placed  at  a measured 
height  above  the  table  so  as  to  have  the  bonded  strip's  free  edge 
perpendicular  to  the  middle  bar's  surfaces.  Use  of  an  angled  square 
aligned  the  three  bars,  adhesive  was  applied  and  the  two  free  edges 
of  the  strips  were  bonded  to  opposite  sides  of  the  middle  bar.  After 
initial  cure,  the  holding  plate,  as  seen  in  Figure  9,  was  fixed  to  the 
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system  to  prevent  movement  and  the  whole  assembly  was  cured  at  1A0°F 
in  a large  walk-in  environmental  room  located  in  Texas  A&M  University's 
McNew  Laboratory. 

Loading  fixtures  were  designed  so  as  to  hold  the  sample  system 
and  be  compatible  with  the  Instron  Machine.  As  in  Figure  11,  the 
samples  could  be  pulled  in  tension  or  shear.  The  tension  was  applied 
through  both  the  samples  at  one  time.  The  samples  were  pulled  in 

simple  shear  in  this  so-called  double  lap  shear  test  (i.e.  the  middle 
bar  stays  fixed  and  the  two  outer  bars  move  up) . In  order  to  perform 
a combination  of  tension  and  shear,  the  sample  system  was  pulled  to 
the  desired  tension  level,  the  holding  plates  were  lightly  screwed  into 
the  outer  bars,  C-clamps  were  applied  to  prevent  rotation  of  the  plates 
upon  tightening,  and  the  plates  were  tightened  with  the  use  of  wrenches. 
Two  holding  plates  were  used,  one  for  each  side,  and  were  designed  to 
hold  the  system  rigidly  in  all  shear  tests,  with  or  without  pretension. 
The  tensioned  sample  system  was  removed  from  the  Instron  and  the  load- 
ing fixture  rotated  90°.  With  the  holding  plates,  the  outer  two  bars 
remain  parallel  to  each  other  when  the  outer  bars  are  pulled  upwards. 
Figure  11  shows  the  described  apparatus  without  the  plates  to  allow  easy 
inspection  for  the  reader. 

Experimental  Considerations 

In  any  experimental  effort  care  must  be  exercised  to  insure  that 
the  intended  experimental  model  and  the  actual  experimental  model  are 
the  same.  Sometimes,  due  to  physical  difficulties,  modifications  to 
the  model  must  be  made  and  then  corrected  for. 
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The  bonding  of  plexiglas  to  stainless  steel  was  no  easy  matter. 

The  two  dissimilar  materials  could  be  bonded  together  with  various 
adhesive  systems,  but  the  high  stress  concentrations  at  the  corners 
of  the  strip  caused  premature  failure.  This  problem  was  eliminated 
by  bonding  0.06  inch  thick  plexiglas  plates  to  the  surfaces  of  the 
stainless  steel  bars  with  an  epoxy  system.  This  epoxy  consisted  of 
40%  Versamid  125  hardener  and  60%  R-815  resin  (by  volume)  as  supplied 
by  the  King  Chemical  Company.  This  system  was  cured  for  24  hours  at 
room  temperature.  The  plexiglas  surfaces  and  edges  were  roughened 
for  bonding  using  grade  200  silicon  paper  and  cleaned  with  acetone, 
M-Prep  Conditioner  A and  B were  used,  then  Freon  TF  was  sprayed  on 
and  wiped  clean  with  tissue.  These  materials  were  obtained  from  Micro- 
measurements. With  these  clean  surfaces  and  edges,  the  biaxial  strips 
were  then  bonded  perpendicular  to  the  plexiglas  plates  with  the  use 
of  Daybond  Thickened  Cement  manufactured  by  Dayton  Plastics  Inc. 

This  allowed  the  bonding  to  be  essentially  cohesive  due  to  the  nature 
of  the  cement  and  to  have  no  excess  that  would  mar  the  orthogonal 
arrangement  of  plate  and  strip.  Then  elevated  temperature  was  used 
for  curing,  as  described  earlier,  for  at  least  three  hours. 

The  elevated  temperature  (at  approximately  140°F)  was  low  enough 
to  prevent  any  permanent  damage  to  the  samples  such  as  the  dimensional 
changes  at  high  temperature  predicted  by  the  manufacturer  of  the  GM 
plexiglas  used.  The  material  at  140°F  was  well  below  the  glass  trans- 
ition temperature  of  plexiglas,  nominally  220°F.  This  temperature  is 
defined  as  the  temperature  above  which  the  physical  mechanisms  of 
deformation  within  the  polymer  start  to  change  due  to  a significant 
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change  in  free  volume [ 1 23 ] . Since  the  tests  were  conducted  at  75°F 
at  50%  relative  humidity  (the  samples  were  at  these  conditions  for  at 
least  twelve  hours,  prior  to  testing)  viscoelastic  behavior  was 
negligible.  Strain  gage  heating  was  neglected  since  others  have  found 
this  to  be  negligible  [10].  Temperature  and  humidity  were  recorded 
in  all  tests  to  check  any  affects,  and  none  were  found.  In  order  to 
check  for  rate  dependence,  some  tests  were  made  at  three  displacement 
rates  (0.002,  0.02,  0.2  in/min).  Little  rate  dependence  could  be 
found  at  these  test  conditions. 

Due  to  the  change  in  the  stress  field  and  boundary  conditions 
of  the  biaxial  strip  because  of  the  addition  of  the  plexiglas  plates, 
an  analysis  was  undertaken.  Viewing  the  strip  and  plates  as  an  I-beam, 
(seen  from  the  end  of  the  strip)  a plane  strain  finite  element  analysis 
was  made.  (The  program  used  will  be  described  in  a later  section). 
Checking  ail  conditions  where  the  different  strips  were  greater, 
equal  to,  and  smaller  in  thickness  than  the  pre-bonded  plates,  we  found 
by  adding  approximately  one  half  of  each  plate  thickness  to  the 
strip's  height,  the  correct  biaxial  stress  field  was  achieved.  Thus 
the  strips'  heights  were  considered  to  be  1.06  inches. 

When  conducting  any  multi-axial  loading  experiment,  proper 
alignment  is  critical.  Alignment  was  kept  as  correct  as  possible  in 
the  bonding  precedure  described  earlier.  This  and  the  loading  align- 
ment were  checked  by  a series  of  tests  with  strain-gaged  sample 
systems.  The  gages  described  earlier  were  fixed  on  both  ends  of  the 


strip,  one  and  one-half  inches  away  from  the  edge  (to  avoid  end 
effects)  on  opposite  sides  of  the  strip.  The  system  had  four  gages. 
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two  to  a strip.  This  can  be  seen  in  Figure  11,  The  gages  were  moni- 
tored with  the  use  of  Vishay/Ellis  - 10  Portable  Strain  Indicators. 

Very  small  amounts  of  bending  occurred  in  loading  in  tension,  which 
was  eliminated  during  the  placing  of  holding  plates.  After  taking 
the  tensioned  sample  system  from  the  Instron  Machine,  strains  were 
still  recorded.  Less  than  5%  decay  in  tensile  strain  was  achieved 
if  the  holding  plates  were  properly  fixed.  Upon  shear  loading, 
negligible  slipping  was  recorded  by  the  strain  indicators. 

Tests  were  made  to  check  how  close  the  man-made  cracks  resembled 
natural  flaws.  By  scribing  a plexiglas  strip  on  one  side  with  a tool 
and  gently  tapping  the  other  side  with  a chisel  and  hammer,  a very 
natural  crack  can  be  made.  Due  to  the  difficulty  in  producing  a number 
of  these  cracks,  only  a few  strips  were  tested  with  these  flaws. 

Little  difference  could  be  detected  between  test  results  with  these 
more  natural  cracks  and  the  machined  crack's  test  results. 

Testing 

Much  of  the  test  procedure  has  been  described  already,  but  a 
short  unified  presentation  will  be  presented  here.  The  strips  were 
measured  with  the  use  of  calipers  and,  as  nearly  as  possible,  matched 
strips  were  used  together.  After  bonding  to  the  steel  bars,  a test 
ratio  of  tension  to  shear  load  was  decided  upon,  and  the  test  conducted. 
Tension,  tension/simple  shear  and  simple  shear  tests  were  conducted 
with  the  center  crack  samples.  With  the  off-center  crack  samples, 
only  tension  and  simple  shear  tests  were  made.  Table  1 is  a summary 
of  this  test  program.  All  the  previously  described  experimental 
considerations  were  taken  into  account. 


Table  1.  Summary  of  Fracture  Tests 
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The  samples  were  pulled  to  failure  and  the  failure  loads  (tension, 
shear,  tension  and  shear)  were  recorded.  Since  not  all  four  crack 
tips  Initiated  at  once,  four  sets  of  loads  resulted.  The  cranes 
started  quite  suddenly,  propagated  quickly  in  a defined  direction,  and 
stopped  at  the  interface  between  sheet  and  strip.  Further  increases 
in  the  applied  displacements  caused  slow  propagation  down  the  inter- 
face. Little  out  of  plane  bending  resulted  as  the  crack  tips  initially 
propagated;  but  after  much  propagation  down  the  interface,  bending 
occured  so  the  tests  were  halted.  Examples  of  the  typical  load  his- 
tories of  the  tests  are  shown  in  Figure  12. 

Data  Reduction 

In  order  to  calculate  the  correct  initial  stress  intensity 
factors  (at  failure),  a plane  stress  finite  element  analysis  was  conducted 
of  the  biaxial  strip.  Modeling  the  strip  (3  x 1.06  x 1 inches) 
was  accomplished  with  cracks  of  lengths  0.25,  0.50,  0.75,  1.25,  1.50, 

2.00,  and  4.00  inches  in  length.  (The  computer  programs  used  will  be 
described  in  a later  section).  The  tension  test  was  modeled  by  displac- 
ing rigidly  the  long  edge  boundaries  perpendicular  to  the  crack  (holding 
displacement  parallel  to  the  crack  fixed).  The  simple  shear  test  was 
modeled  by  displacing  rigidly  the  long  edge  boundaries  parallel  to 
the  crack  (holding  displacement  perpendicular  to  the  crack  fixed). 

These  two  tests  were  modeled  over  the  range  of  c rack  lengths  used. 

The  load  for  a unit  displacement  was  found  in  all  the  cases  by  inte- 
grating the  stresses  over  the  boundary.  The  analysis  also  confirmed 
that  the  strain  gages  used  for  checking  alignment  were  correctly 
placed  to  avoid  end  effects.  The  stress  intensity  factors  were  plotted 


Figure  12.  Examples  of  loading  histories  (resultant  load  vs  cross-head  displacement) 
used  in  fracture  tests. 
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versus  crack  length/strip  height  as  in  Figures  13  and  14.  The 
solid  and  dashed  curves  represent  equations  for  a crack  in  an  all- 
around  infinite  sheet  (solid  line)  [36],  and  a crack  in  an  infinitely 
long  strip  [121],  The  finite  element  results  show  the  close  agreement 
with  theory.  The  slightly  lower  values  for  the  stress  intensity 
factors  are  believed  due  to  the  strip  not  being  infinitely  long.  The 
finite  element  computed  load  versus  crack  length/strip  height  is 
plotted  as  in  Figures  15  and  16.  Then  the  ratio  of  stress  intensity 
factor  to  load  was  plotted  versus  crack  length/strip  height  as  in 
Figures  17  and  18.  Using  these  last  figures,  knowing  the  crack  length 
of  the  sample,  and  the  loads  causing  crack  growth,  the  correct  stress 
intensity  factor  could  be  calculated  for  any  amount  of  grip  displacement. 

The  angles  of  crack  propagation  for  combined  loadings  were 
initially  measured  directly  from  the  samples.  As  this  was  difficult, 
later  samples  were  removed  from  the  bars  by  a chisel,  laid  upon  paper 
Vand  crack  path  traced.  The  initial  angle  of  propagation  was  used  for 
comparison  with  theory.  It  was  found  that  the  crack  propagated  along 
a straight  line  for  a significant  distance,  then  curved  to  meet  the 
boundary  as  shown  in  the  section  on  f ractography . This  distance  decreased 
for  increased  levels  of  tension  in  combined  mode  testing.  The  angle 
between  this  initially  straight  line  and  the  crack  was  used  as  the 
initial  propagation  angle  . 

Results 

The  results  were  plotted  in  such  a way  that  they  could  be  compared 
easily  with  the  current  theories  describing  crack  trajectories. 

Figure  19  shows  the  ratio  of  critical  stress  intensity  factors  versus 


Figure  14.  Sliding  mode  stress  intensity  factor  vs  crack  length/ 
strip  height  for  unit  displacement. 


Figure  15.  Tension  load  vs  crack  length/strip  height 
for  unit  displacement. 


Figure  19.  Crack  propagation  angle  vs  stress  intensit 
theoretical  maximum  energy  release  rate  [9 
data  points. 
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crack  propagation  angle.  These  points  were  compared  to  the  various 
theories.  The  theoretical  energy  release  rate  results  of  Knauss  and 
Palaniswamy  [92]  compared  quite  well  in  certain  portions  of  the  data 
• and  proved  to  be  best  overall.  (Referring  to  Figure  3,  it  is  found  that 

the  maximum  stress  criterion  based  on  the  singular  solution  predicts 
values  which  fall  along  the  lower  edge  of  the  scatter  band  in  Figure  19.) 
The  dashed  line  shows  the  area  left  incomplete  by  the  theoretical 
analysis,  except  for  the  pure  mode  II  problem. 

A plot  of  the  relationship  between  the  stress  intensity  factors  is 
shown  in  Figure  20.  These  points  show  the  problem  inherent  in  most 
fracture  data;  viz.,  a fair  amount  of  scatter.  The  curve  was  generated 
| I from  predictions  made  by  the  energy  release  rate  criterion.  The  criterion 

fits  the  data  reasonably  well  and  the  plot  serves  as  a failure  curve  for 
plane  problems.  However,  the  maximum  stress  criterion  is  slightly  better 
at  small  angles  (cf.  Figure  6). 

Referring  to  Figure  19,  it  is  seen  that  the  angle  of  propagation 
should  be  predicted  in  problems  with  low  K^/K^  in  order  to  complete  the 
solid  curve;  this  will  be  investigated  by  the  use  of  energy  release  tech- 
niques with  finite  element  codes  in  a later  section.  The  reason  why  a 


greater  load  is  needed  to  cause  propagation  down  the  interface  will  also 
be  examined.  Also,  the  ability  to  predict  the  angle  change  as  the  crack 
approaches  the  interface  will  be  studied. 


Sliding  vs  opening  stress  intensity  factors  for  theoretical 
maximum  energy  release  rate  [92]  and  experimental  data  points 


Various  interesting  observations  can  be  made  concerning  the 
I ■ crack's  path.  The  tests  with  the  center  cracks  in  either  tension  or 


shear  without  tension  provided  little  new  information.  The  cracks 
either  propagated  fairly  straight  down  the  center  or  straight  to  the 


boundary.  The  off-center  cracks  in  tension  and  shear  acted  similarly. 
The  off-center  cracks  in  tension  propagated  along  a slight  curve  to 
the  center  of  the  strip  and  then  down  the  middle  to  the  ends  of  the 
strip  as  predicted  by  Erdogan  [18].  As  described  earlier,  the  center 
cracks  under  combined  tension  (K^)  and  shear  (K^)  propagated  on  a 
straight  line  at  a well-defined  angle  until  approaching  the  interface. 
Depending  on  the  percentage  of  tension,  the  crack  would  curve  to  meet 
the  interface.  Pictures  taken  with  a bel lows-camera  system  shows  this 
phenomenon  in  Figure  21. 

Using  an  optical  microscope,  a Leitz  instrument  with  magnification 
up  to  220X,  the  angle  of  propagation  next  to  the  intitial  tip  was 
examined.  The  surfaces  of  the  cracks  were  also  examined.  The  shallow 
depth  of  field  of  this  equipment  made  examination  difficult. 

A more  sophisticated  instrument,  the  JEOL  JSM-U3  scanning  electron 

L , microscope  (SEM)  was  also  used.  This  SEM  has  a 100  X resolution  with 

r 

[*  magnification  from  22X  to  100,000X.  Observations  could  be  made  of  the 

samples'  surfaces  and  recorded  by  the  use  of  photographv.  The  SEM 

was  operated  in  the  secondary  electron  image  mode,  with  various  control 

settings  and  magnifications.  This  mode,  where  the  electrons  detected 

are  emitted  from  the  sample  or  its  conductive  gold-palladium  coating, 

I | 

possesses  a great  ability  for  detection  of  details  with  a great  depth 

■ i 

A 
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of  field  [124].  This  instrument  was  used  to  investigate  the  cracks' 
angles  and  surfaces  as  had  others  in  investigating  composite's  surfaces 
[125,126]. 

In  examining  the  initial  angle  up  to  magnifications  up  to  10,000X 
it  was  found  that  there  was  a very  small  region  where  the  crack  tip 
was  not  perpendicular  to  the  surface  of  the  strip.  This  was  roughly 
proportional  to  the  amount  of  curvature  in  the  initial  crack  front, 
prior  to  loading.  If  the  initial  crack  front  was  curved,  the  propa- 
gating front  grew  so  as  to  align  itself  perpendicular  to  the  sheet 
surface;  the  more  initial  curvature,  the  greater  the  distance  for 
alignment.  This  region  of  alignment  was  very  small,  usually  less 
than  the  sheet  thickness.  This  is  shown  in  Figure  22.  The  angle  of 
propagation  used  was  the  angle  between  the  initial  crack  and  the  crack 
path  after  alignment.  After  examining  a number  of  crack  tips,  no  correla- 
tion could  be  made  between  the  very  small  curvature  at  the  initial 
crack  tip  shown  in  Figure  22  and  the  craze  zone  size  reported  for 
plexiglas  [121].  No  evidence  of  ductility  or  crazing  marring  the  crack 
surface  could  be  detected. 

Investigating  the  crack  surfaces  revealed  the  same  appearance 
reported  as  others  [127,128].  At  slow  crack  speeds,  a large  number 
of  parabolic  curves  roughen  the  surface,  as  speed  increases  the  surface 
is  smooth,  and  at  speeds  approaching  the  limiting  elastic  wave  velocity 
only  lines  parallel  to  crack  direction  can  be  seen.  At  the  initial 
crack  tip  after  the  crack  front  is  aligned  the  surfaces  showed  the 
slow  crack  velocity  markings.  The  markings  then  decreased  until  onlv 
the  grooves  remained,  showing  increased  speed  as  the  crack  propagated 
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to  the  interface.  This  can  be  seen  in  Figure  23.  Near  the  poin 
where  the  crack  trajectory  curved,  the  markings  again  appeared  iridicat 
ing  a slowing  process.  More  and  more  markings  appeared,  indicating 
crack  arrest  as  the  crack  approached  an  interface,  as  shown  in  Figure 

24. 

This  examination  of  the  crack' s surface  showed  that  a properly 
made  crack  is  essential  in  fracture  tests.  Investigation  of  surface 
features  also  provided  evidence  of  crack  arrest  as  the  crack  tip 


approached  a bonded  edge. 
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FINITE  ELEMENT  ANALYSIS 
General 

Certain  numerical  analyses  were  undertaken  to  supplement  existing 
analytical  results  and  investigate  problems  encountered  in  the  experi- 
mental effort.  Due  to  its  flexibility  in  solving  complicated  problems, 
the  finite  element  method  was  employed.  Two  different  programs  v ere 
used.  After  reviewing  the  finite  element-fracture  techniques  and  using 
some  of  them,  it  was  determined  that  codes  containing  special  cr;  ck 
elements  would  be  the  best  for  the  mixed  mode  plane  problem.  It  was 
further  decided  that  the  hybrid  elements  based  on  Pian's  work  [62] 
would  be  the  best  for  our  purposes  due  to  the  reasons  described  in  a 
previous  section. 

The  first  program  used  was  a modified  version  of  the  code  listed 
in  Desai  and  Abel's  book  [129].  This  is  a plane  strain/stress  computer 
code  which  is  limited  to  linear,  elastic  isotropic  bodies.  The  elements 
used  are  quadrilaterals  composed  of  four  linear  displacement  triangles 
(with  condensation  of  the  internal  node)  and  linear  displacement 
triangles.  Modifications  were  two-fold.  First,  the  capacity  of  the 
program  was  extended  by  the  use  of  larger  arrays  and  more  scratch 
files,  to  allow  larger  problems  to  be  solved  with  more  elements. 
Secondly,  the  fracture  analysis  capability  was  enhanced  by  calculating 
and  printing  the  strain  energies  of  the  elements  and  body,  and  by 
the  addition  of  a crack  element  to  the  element  library.  This  ele- 
ment is  the  hybrid  crack  element  presented  in  the  literature  [62- 
64,  130],  It  is  compatible  with  the  other  elements  of  the  code. 
Extensive  check  problems  were  investigated  to  learn  the  capabilities 
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and  limits  to  this  modified  program  called  CRACK. 

In  order  to  solve  larger  problems  more  accurately,  the  Texas 
Grain  Analysis  Program,  TEXGAP  [131],  was  used,  which  is  a finite 
element  program  used  for  analyzing  solid  rocket  propellant.  The  ver- 
sion of  the  code  employed  here  is  for  linearly  elastic  plane  or  axisym- 
metric  bodies.  The  elements  used  were  reformulated  isotropic  quad- 
rilaterals composed  of  four  quadratic  di splacement,  linear  H triangles 
(allowing  Poisson's  ratio  to  be  0.5  for  plane  problems),  non-reformu- 
lated  quadratic  displacement  subparametri c quadrilaterals,  reformulated 
isotropic  quadratic  displacement, linear  H triangles,  and  hybrid  crack 
elements . 

In  support  of  the  experimental  effort,  TEXGAP  was  used  to  analyze 
the  affect  of  the  bonding  geometry  on  the  biaxial  stress  field.  CRACK 
was  used  for  the  determination  of  the  stress  intensity  factors  with 
some  checks  made  with  TEXGAP.  Both  programs  were  well-suited  for  the 
analysis  of  fractured  bodies. 

Biaxial  Strip  Ai^  lysis 

A study  was  made  to  determine  if  the  experimental  crack  growth 
behavior  could  be  predicted.  This  behavior  was  divided  into  three 
phases.  The  original  configuration  of  a centered  crack  was  used  to  deter- 
mine the  initial  crack  propagation  angle.  This  analysis  was  then  used 
to  predict  the  change  in  the  angle  of  propagation  as  the  crack  grew. 

The  next  phase  was  to  investigate  why,  after  the  crack  arrested  at 
the  interface.  Increased  applied  displacement  was  necessary  to  cause 
further  propagation  along  the  Interface.  These  three  phases  were 
analyzed  with  the  use  of  the  finite  element  codes  mentioned  previously. 
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Due  to  the  limitations  pointed  out  earlier  in  the  crack  propaga- 
tion theories,  an  examination  of  the  maximum  energy  release  criterion 
with  the  use  of  finite  elements  was  also  undertaken.  The  experimental 
biaxial  strip  was  modeled  with  the  CRACK  program.  Two  types  of  analysis 
were  made,  one  with  the  use  of  crack  elements  and  one  with  only  conven- 
tional elements. 

With  the  use  of  two  crack  elements  (one  for  each  crack  tip)  for 
a biaxial  strip  configuration,  much  experience  had  been  gained  in  the 
analyses  conducted  in  support  of  the  experimental  work.  Additional 
test  problems  were  conducted  to  check  on  the  degree  of  mesh  refinement 
needed  to  calculate  total  strain  energy.  For  a uniaxial  bar  problem, 
increasing  the  number  of  elements  from  five  to  forty  brought  the  cal- 
culation of  strain  energy  from  99.4%  to  99.9%  of  the  exact  value  of 
strain  energy.  Since  a comparison  of  strain  energies  for  a constant 
mesh  was  to  be  made,  this  error  was  felt  to  be  negligible.  To  check 
the  affect  of  the  crack  elements  on  the  strain  energy,  two  mesh 
geometries  were  used.  The  refined  model  used  crack  elements  of  one- 
half  the  size  of  the  other  mesh.  Twice  as  many  conventional  elements 
were  used  about  the  crack  elements  in  the  refined  mesh,  as  compared 
to  the  other  mesh.  Negligible  difference  was  obtained  in  the  angle 
of  the  greatest  energy  release. 

Using  finite  element  models  of  strips  of  the  same  proportions  as 
the  experimental  ones,  the  strain  energies  of  center  cracked  and  incre- 
mently  cracked  biaxial  strips  were  analyzed.  Varying  ratios  of  applied 
normal  and  shearing  displacements  (corresponding  to  different  ratios 
of  stress  intensity  factors)  along  the  top  and  bottom  of  the  two  strip 
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models  were  applied  through  a rigid  boundary.  In  the  second  geometry, 
the  crack  increment's  length  (0.01  < AA  < 0.1)  and  angle  were  varied  to 
search  for  the  maximum  energy  release  as  calculated  by: 


G = 


W initial  - W incremented 


2 Aa 


(16) 


When  G was  maximized  for  a particular  angle,  this  was  assumed  to  be 
the  angle  of  propagation. 

In  using  this  procedure,  it  was  found  that  only  a certain  small 
range  of  increment  lengths  would  produce  a maximum  for  a particular 
displacement  ratio.  The  angle  of  propagation  was  constant  for  this 
range . As  the  crack  propagation  angle  increased  (corresponding  to  lower 


K^/K.^  ratios)  the  maximum  was  harder  to  determine.  It  was  felt  that 


this  was  due  to  the  crack  element's  inability  to  handle  crack  tips 
that  were  too  far  from  the  element's  center.  Off-setting  the  element 
as  in  Figure  25,  helped  correct  this  a little. 

Due  to  these  problems,  another  mesh  configuration  using  only 
conventional  elements  was  used  as  in  Figure  25.  This  mesh  was  used 
to  check  and  correct  some  of  the  previous  answers,  especially  for  cases 
involving  low  K ^ / K ^ ratios.  In  this  geometry,  the  crack  tip's  nodal 
point  was  changed  for  various  angles  and  increment  lengths  to  determine 
the  maximum  value  of  G.  This  method  was  less  sensitive  to  changes  in 
increment  lengths  and  allowed  a more  reliable  method  of  determining 
the  maximum  energy  release,  G,  and  the  angle  of  propagation,  9. 

The  question  of  using  a finite  increment  for  calculating  G must 
be  faced.  Due  to  the  finite  element  representation,  some  finite  value 
has  to  be  used.  The  ratios  of  maximum  G for  combined  mode  to  the  G 
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for  opening  mode  for  this  analysis  were  calculated.  Comparison  to  the 
values  from  Hussian's  et  al . [95]  theoretical  analysis  for  a virtual 
(approaching  zero)  increment,  the  ratios  were  in  close  agreement.  The 
finite  element  results  were  therefore  below  those  given  in  [92]  (cf. 

Figure  6) . 

The  results  of  examining  this  first  phase  of  crack  propagation  are 
presented  in  Figure  2b.  The  initial  angle  of  propagation  is  plotted 
against  the  ratio  of  stress  intensity  factors  for  the  original  crack. 

This  curve,  while  calculated  from  data  for  a specific  geometry, 
can  be  used  for  any  two-dimensional  crack  problem.  It  can  also  be 
compared  to  the  experimental  results  of  this  report  as  in  Figure  19, 
as  well  as  the  results  of  other  theories.  The  results  for  this  finite 
element  analysis  are  slightly  lower  than  that  of  Knauss  and  Palaniswamy ' s 
analytical  results  [92],  but  agree  well  with  Coughlan  and  Barr's  [96] 
finite  element  results.  The  results  also  provide  9 over  the  range  of 
values  of  K|/K,  not  fully  covered  in  other  researcher's  analyses.  The  high 
range  of  values  of  K^/K0  were  not  investigated  due  to  the  agreement 
already  found  in  the  literature. 

With  the  first  phase  of  crack  propagation  completed,  the  analysis 
of  the  crack  that  has  propagated  part-way  to  the  clamped  boundary  was 
made.  Two  cases  were  examined,  based  on  the  initial  angle  of  propa- 
gation. These  angles  were  75"  and  55°  corresponding  to  initial  stress 
intensity  factor  ratios,  K^/K^,  of  0.  and  1.1.  The  program  TEXGAP 
was  used  in  this  analysis.  The  crack  geometry  used  for  the  finite 
length  angled  crack  tip  is  shown  in  Figure  27.  The  other  crack  tip 
was  also  incremented  but  not  modeled  with  a crack  element.  The  crack 


igure  26.  Crack  propagation  angle  vs  stress  intensity  factor  ratio  for  finite 
element  results  compared  against  other  criteria. 
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element  size  was  varied  as  also  the  number  of  elements  behind  anti  ahead 
of  the  crack  element  were  increased.  Little  change  in  answers  were 
found.  The  displacements  required  to  produce  the  initial  stress  inten- 
sity factors  were  also  applied  to  the  biaxial  strip  with  the  angle 
crack.  The  new  stress  intensity  factors  for  this  geometry  were  then 
obtained. 

These  new  stress  intensity  factors  formed  a ratio  which 

could  then  be  used  to  predict  the  direction  of  crack  propagation.  The 
predicted  initial  angle  of  propagation  (0=55°)  produced  a new  ratio  of 
2.9  with  a new  propagation  angle,  $=32°,  to  the  angled  crack  line  of 
propagation;  the  angles  shown  in  Figure  27b  are  52°  and  24°,  respectively. 
The  mode  II  problem  (9=75°)  had  a new  ratio  of  18.5  with  a new  propagation 
angle  of  <£=5°;  from  Figure  27a,  these  angles  are  70°  and  6°,  respectively. 
These  new  propagation  paths  did  not  align  themselves  perpendicular  to  the 
direction  of  global  maximum  principal  stress,  indicating  the  need  of  a 
local  criterion  to  govern  crack  direction. 

The  next  phase  of  the  crack's  behavior  to  be  analyzed  was  the 
matter  of  arrest  and  propagation  at  the  interface.  Two  cases  were 
considered:  cracks  with  propagation  angles  of  65°  (combined  mode) 

and  75°  (mode  II).  The  program  TEXGAP  was  used  with  no  crack  elements. 

The  total  energy  release  rates  at  the  initial  center  crack  were  cal- 
culated for  a certain  increment  length.  The  total  energy  release 
rates  were  in  approximate  agreement  with  Hussain's  et  al ■ [95]  energy 
release  rate  mentioned  earlier.  The  crack  was  then  allowed  to 
propagate  to  the  interface  by  means  of  disconnecting  elements.  Ihe 
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energy  of  this  configuration  was  determined  for  the  same  applied  displace- 
ments as  the  initial  crack.  The  energy  release  rate  along  the  interfaces 
was  then  calculated  for  the  same  increment  length  as  for  the  initial 
crack  to  avoid  any  effect  of  numerical  errors. 

For  specified  applied  displacements,  the  total  energy  release  for 
the  center  crack  was  2.1  times  the  energy  release  for  propagation  along 
the  interface  for  the  combined  mode  case.  The  sliding  mode  case  had  a 
ratio  of  3.4.  Assuming  the  initial  crack  and  interface  crack  propagates 
at  the  same  critical  value  of  the  total  energy  release  rate,  the  crack 
arrests  at  the  interface. 

In  order  to  calculate  the  displacement  required  to  propagate  a crack 
along  the  interface,  it  is  helpful  to  draw  upon  the  linearity  of  the 
problem.  Namely,  in  view  of  this  linearity,  the  relation  between  the 
total  energy  release  rate  and  the  applied  displacements  can  be  written: 


2 2 
G = y u + y u v + y v 
11  o 12  o o 22  o 


where  Y,-,,  Yim  and  are  independent  of  u and  v . In  the  case  of  the 
11  1 z zz  o o 

sliding  mode,  the  ratio  of  displacement  needed  to  propagate  the  crack 
along  the  interface  to  that  for  the  center  crack  is  easily  found  to  be 
/3.4  = 1.84  since  v^  = 0 in  this  case.  Table  2 shows  the  comparison  of 
the  analytical  and  experimental  values  of  increased  displacement  needed 
to  cause  further  propagation. 


Table  2.  Percent  of  Initial  Displacement  Needed  to  Cause  Interface 
Crack  Propagation. 


Analytical 

Experimental 

Mixed  Mode 

188% 

175-200% 

Shearing  Mode 

184% 

175-200% 

The  range  of  values  for  the  experimental  results  is  due  to  the  varia- 
bility of  the  cohesive  bond  between  the  plexiglas  plates  and  sheets 
as  described  in  the  experimental  section  of  this  study. 

From  the  previous  discussion,  it  is  believed  that  the  crack  path 
can  be  predicted  in  the  biaxial  strip  model  of  a composite.  For  this 
model  with  very  stiff  or  rigid  fibers,  it  has  been  demonstrated  that 
cracks  will  propagate  toward  the  boundary  provided  some  of  the  growth 
is  of  the  second  mode  (sliding  mode).  It  also  has  been  shown  that  this 
propagation  path  can  be  predicted.  The  crack  arrests  at  the  boundary 
and  increases  in  applied  displacements  must  be  made  to  cause  further 
crack  growth.  With  a crack  arrest  mechanism  described  for  composites 
with  combined  mode  and  shearing  mode  microcrack  propagation,  a crack 
arrest  mechanism  for  opening  mode  microcrack  propagation  must  be 
explored.  This  will  be  done  by  considering  fiber  flexibility. 


Simulated  Composite  Analysis 

A natural  extension  of  the  previous  analysis  is  to  go  to  a model 
with  flexible  fibers.  The  motivation  for  this  is  to  calculate  the 
energy  release  rate  due  to  fiber  bending  and  to  show  the  affect  on 
crack  growth,  especially  on  cracks  in  an  opening  mode  situation.  In 


this  simple  model  of  a composite  subject  to  loads  at  infinity,  the  fibers 
just  above  and  below  the  crack  tend  to  respond  as  beams  on  an  elastic 
foundation.  With  no  crack  present,  the  fibers  would  just  separate 
slightly.  With  a crack,  the  crack  is  equivalent  to  an  internal  pressure 
(equal  and  opposite  to  the  stress  in  an  uncracked  body)  acting  between  the 
fibers  of  an  externally  unloaded  specimen,  thus  causing  local  bending. 

This  approach  was  taken  in  a previous  analysis  using  a simple 
model  of  a beam  on  an  elastic  foundation  [8].  For  the  case  of  a uni- 
formly loaded  fixed-ended  beam  on  an  elastic  foundation  [132],  the 
deflection  and  loading  was  used  to  calculate  the  strain  energy  release 
rate  due  to  beam  bending.  The  relationships  between  energy  release 
due  to  beam  bending,  crack  length,  and  fiber  rigiditv  were  explored. 

A simple  but  more  realistic  model  of  a composite  exhibiting  fiber 
bending  was  needed.  A finite  element  representation  of  the  geometry 
in  Figure  28  was  examined  using  TEXGAP.  The  cases  were  examined  either 
as  pure  opening  mode  or  as  pure  sliding  mode  where  uniform  displace- 
ments were  applied  at  the  fixed  boundary.  Only  the  right  half  of  the 
problem  was  modeled  with  the  use  of  the  crack  element,  taking  advantage 
of  the  symmetry  of  the  problem.  Test  cases  of  different  dimensions 
were  run  to  ensure  the  model  was  long  enough  to  be  infinitely  long. 

In  this  case,  a size  of  5.  x 40.  was  needed  for  the  total  model. 

The  problem  was  analyzed  as  if  it  were  in  plane  stress;  for 
plane  strain  (plates  in  bending)  the  results  would  be  qualitatively 
the  same  and  plane  strain  results  can  be  extracted  from  the  plane  stress 
analysis  by  an  appropriate  change  in  the  material  constants.  The 
material  constants  used  in  this  analysis  are  those  typical  of  a graphite 


.35,  hnj  = 1.).  Actual  size  of  model  is  5.  x 40. 
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fiber-reinforced  epoxy  matrix.  For  equal  fiber,  h^,  and  matrix,  h , 
heights  in  an  idealized  circular  fiber  composite  [81,  the  fiber  volume 
fraction,  V^,  is  computed  by 


which  gives  a value  of  0.71,  not  unreasonable  for  advanced  composites. 

The  deflection  in  the  fiber's  y-direction  about  the  crack  tip  for 
a uniform  vertical  grip  displacement  v^  are  presented  in  Figure  29. 

As  can  be  expected,  no  fiber  bending  occurred  in  the  shear  case  and 


it  will  not  be  presented.  For  small  cracks,  the  center  line  deflec- 
tions of  the  fiber  are  not  large,  but  their  bending  extends  almost  ten 


times  the  crack  length.  With  increases  in  the  initial  crack  length,  or  as 
the  crack  grows,  the  deflection  gets  much  larger,  but  the  length  of  the 
fiber  affected  only  slowly  increases.  As  the  crack  approaches  a length 
of  ten  times  the  fiber  height,  the  portion  of  the  fiber  above  the  crack 
deflects  nearly  as  much  as  the  applied  displacement.  (This  is  due  to 
the  use  of  a uniform  grip  displacement.)  The  effect  along  the 
fiber  is  only  over  twice  the  crack  length.  As  the  crack  gets  larger, 
the  deflection  does  not  get  larger,  but  the  length  of  the  fiber  affected 
beings  to  slowly  increase. 

Values  of  the  total  energy  release  rate  G . were  calculated 
from  the  stress  intensity  factors  (from  the  finite  element  analysis) 
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Figure  30.  Energy  release  rates  vs  y for  opening  mode.  ( -» .finite  element  solution' 
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where  k is  Che  foundation  stiffness  and  I is  the  fiber's  area  moment  of 
inertia.  The  value  of  k for  this  analysis  and  for  an  actual  composite 
could  be  treated  as  a free  parameter  used  to  correlate  theory  and  experi- 
ment. For  the  circular-fiber  composite,  an  estimate  was  made  for  k 
in  [8], 


E h, 
m f 


E /V, 

m f 


(1+v  ) h (1+v  ) l-/v7  K J 

mm  m f 

A better  estimate  for  k can  be  calculated  with  the  use  of  an  equation 

•k 

derived  in  the  previously  referenced  research  [8], 

o 2 he  ,2  2 

C = —P—L  r^ycosh  Y ~ cos  y)  , , 

BI  k (sinh  y + siny) ^ (24) 


Figure  30  shows  how  the  total  energy  release  rate  in  Mode  I is 
strongly  affected  by  fiber  bending.  For  very  small  cracks,  the  results 
are  bounded  by  the  release  rate  of  a crack  in  an  infinite  sheet  given 


As  the  crack  length  increases,  total  energy  release  rate  follows  the  same 
qualitative  shape  as  for  the  energy  release  rate  due  to  fiber  bending 
calculated  in  [8].  For  sufficiently  long  cracks,  the  total  energy 
release  rate  becomes  independent  of  crack  length  as  does  the  bending 
energy  release  rate.  The  magnitudes  of  the  bending  energy  release 
rates  at  these  crack  lengths  could  be  predicted  by  the  equations 


*The  factor  of  2 in  the  brackets  was  inadvertently  omitted  in  the  final 
typing  of  (8). 
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derived  in  [8]  with  a reasonable  value  of  k (490,000  psi,  which  is 
coincidentally  close  to  the  value  of  the  matrix  modulus). 

Figure  31  shows  how  the  total  energy  release  rate  in  the  Mode  II 
problem  increases  without  bound;  there  is  little  effect  of  fiber 
bending.  It  too  is  bounded  for  very  small  cracks  by  the  energy  re- 
lease rate  for  a crack  in  an  infinite  sheet  given  by 
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The  energy  release  rates  for  an  orthotropic  body  were  calculated 
using  the  effective  properties  of  the  material  used  in  the  finite 
element  analysis  and  are  shown  in  Figures  30  and  31.  The  details  of 
this  analysis  can  be  found  in  the  Appendix.  This  was  undertaken  to 
examine  the  relationship  between  the  overall  composite's  and  fiber 
model's  energy  release  rates.  The  equations  for  the  plane  stress 
energy  release  rates  [133]  are 


1 1 an  2 an 


c _ k-  2 a a 1/2  2a  + a 1/2 

Gn  ‘ Si  —it-  Kjii)  +-  ‘ja  ■■“] 

(2)  ' 11  2 11 


(27) 


where  the  a..'s  are  material  constants  for  an  orthotropic  body. 
Using  the  stress  intensity  factors  for  the  crack  in  an  infinite  sheet, 
these  release  rates  predicted  much  larger  values  than  that  of  the  fiber 


bending  model.  This  indicates  the  importance  of  fiber  bending,  and 
that  many  fibers  may  be  affected  by  even  one  crack. 
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Further  Results 

Using  the  results  from  the  rigid  fiber  models,  some  observations 
can  be  made  as  to  crack  growth  and  arrest.  Assuming  that  the  energy 

r* 

per  unit  area,  , required  to  produce  new  surfaces,  is  independent 
of  Loading,  crack  growth  will  occur  when  G^  > 2^  and  arrest  when  G < 
2''".  Considering  the  entire  crack  front,  the  point  of  initiation  of 
crack  growth  will  be  where  this  criterion  is  met. 

In  order  to  determine  how  fiber  bending  affects  crack  extension 
of  an  initially  centered  crack,  was  plotted  vs  y,  in  Figure  32, 

for  equal  fir-field  tension  and  shear  stresses.  The  rigid  fiber  model 
results  are  also  plotted  for  comparison.  For  very  small  cracks,  the 
crack  behaves  as  if  it  were  in  an  infinite  sheet,  with  K^/K^  = 1, 
corresponding  to  a crack  propagation  angle  of  approximately  55°  (from 
Figure  26).  For  somewhat  larger  initial  cracks,  their  behavior  becomes 
affected  by  fiber  bending,  in  that  K^/K^-1.6,  corresponding  to  a 
crack  angle  of  approximately  47°.  The  rigid  grip  solution  with  K^/K^ 
0.37,  indicates  a crack  propagation  angle  of  approximately  65°.  This 
variation  in  angle  will  not  unduly  affect  crack  propagation;  the 
crack  will  still  propagate  to  the  fibers  at  a fairly  steep  angle  and 
will  then  arrest  at  the  fibers.  These  changes  in  stress  intensity 
factor  ratios  due  to  fiber  bending  will  also  have  an  affect  on  crack 
initiation.  The  change  in  ratio  changes  the  point  of  failure  on  the 
^I^II  *nteract*on  curve  as  in  Figure  6.  ’’igher  ratio  would  involve 
greater  applied  stresses  to  cause  crack  propagation. 

Of  great  interest  is  arrest  mechanisms  for  opening  mode  cracks. 

One  possible  mechanism  is  the  propagation  of  cracks  in  adjacent  matrix 
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areas.  The  shear  stresses  produced  in  this  neighboring  matrix,  caused 
by  fiber  bending,  helps  to  cause  cracks  to  grow  from  pre-existing  flaws, 
cracks  of  sufficient  length  to  relieve  the  stresses  at  the  tip  of  the 
initial  crack.  As  seen  in  Figure  33,  the  shear  stress  directly  above 
the  crack  tip  increases  and  then  levels  off  as  the  crack  length  increases. 
The  relative  magnitude  of  the  shear  stress  to  applied  stress  is  small, 
but  a small  amount  of  increase  in  a stress  concentration  can  be  signifi- 
cant in  determining  crack  initiation;  three  dimensional  effects  will 
probably  increase  this  stress  concentration.  Related  crack  skipping  has 
been  observed  experimentally  [8,23].  The  shear  stress  caused  by  fiber 
bending  within  the  matrix  about  the  crack  is  of  some  significance. 

The  shear  stresses  above  the  crack  tip  and  below  the  fiber  are  plotted 
in  Figure  33.  For  an  off-centered  crack,  a shearing  mode  contribution 
can  be  added  to  what  had  previously  been  only  opening  mode  propagation. 
This  will  then  turn  the  crack  toward  the  fibers,  causing  crack  arrest. 

As  the  crack  continues  to  grow,  either  in  the  interface  between 
the  fiber  and  matrix  or  in  the  matrix,  the  growth  will  continue  to 
be  controlled  by  the  neighboring  fibers.  Comparison  of  the  energy 
release  rates  of  the  fiber  model  and  those  for  the  homogeneous  orthotropic 
model,  in  Figures  30  and  31,  shows  that  a cracked  composite  will  re- 
quire much  higher  applied  stresses  for  small  cracks  to  propagate  than 
had  been  predicted  by  using  the  latter  theory.  This  indicates  that  the 
crack  must  be  very  large,  compared  to  the  i her  spacing  and  diameter 
discussed  here,  in  order  to  treat  the  material  as  homogeneous. 

Comparison  of  energy  release  rates  and  stress  intensity  factors 
have  lead  to  some  conclusions  on  mlcrocrack  growth  and  arrest  in 
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fibrous  composites.  The  simple  models  used  here  have  neglected  many 
aspects  of  real  composites  such  as  fiber  randomness,  flaw  size  distri- 
bution and  interface  inhomogeneity.  The  complex,  real  three-dimensional 
problem  was  also  approximated  with  two-dimensional  models.  Even  with 
these  approximations,  it  is  believed  that  these  models  portray  some 
of  the  features  of  microcracking  found  in  the  actual  composite. 
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CONCLUSIONS  AND  RECOMMENDATIONS 

After  examining  the  literature,  it  was  concluded  that  of  all 
existing  criteria,  the  maximum  energy  release  rate  criterion  had  the 
strongest  physical  basis  and  best  correlation  with  data,  especially 
for  brittle  materials.  The  usefulness  of  the  finite  element  method 
in  fracture  mechanics  analysis  has  also  been  demonstrated  in  the 
literature  surveyed.  The  maximum  energy  release  rate  criterion  was 
extended  to  those  combined  mode  fracture  problems  which  are  largely 
composed  of  the  shearing  mode.  This  was  accomplished  with  the  use  of 
the  finite  element  method  and  supported  with  experimental  evidence 
provided  with  a test  geometry  not  previously  used  for  combined  mode 
problems . 

This  completed  criterion  was  then  used  to  predict  the  crack  paths 
and  stress  intensity  factors  interaction  in  long  rectangular  strips 
of  plexiglas  clamped  on  the  long  edges.  The  predictions  compared 
favorably  with  experimental  evidence.  By  correlating  energy  release 
rates  of  cracks  initially  in  the  center  of  the  strip  with  cracks 
already  propagated  to  the  long  edges,  it  was  calculated  and  experi- 
mentally verified  that  significantly  increased  displacements  were 
necessary  to  continue  crack  growth. 

Building  upon  the  biaxial  strip  model  (or  rigid  fiber  model), 
the  affect  of  fiber  flexibility  on  crack  arrest  was  explored.  Fiber 
bending  decreased  the  tendency  of  the  crack  to  propagate  toward  the 
fibers,  but  supplementary  crack  mechanisms  resulted.  These  came  about 
due  to  the  shear  stresses  induced  by  bending;  shear  stresses  in  adjacent 
matrix  regions  possibly  causing  crack  skipping  and  shear  stresses 


in  the  matrix  about  the  crack  causing  propagation  toward  the  fibers, 
both  resulting  in  increasing  the  tendency  for  crack  arrest.  Drawing 
upon  both  the  rigid  fiber  and  flexible  fiber  models  used,  the  nature 
of  crack  growth  from  a small  defect  to  a crack  whose  dimensions  are 
comparable  with  fiber  spacing  was  discussed. 

Examination  of  fractured  surfaces,  obtained  from  the  experimental 
tests,  showed  that  the  direction  of  initial  crack  growth  was  contin- 
uously changing,  reaching  a constant  growth  direction  within  a distance 

which  was  less  than  sheet  thickness  but  much  larger  than  the  estimated 

o 

craze  zone  size  (approximately  6800  A).  The  crack  surfaces  showed 
that  the  crack  velocity  increased  as  the  crack  initially  propagated 
and  then  crack  velocity  decreased  as  the  crack  approached  the  long  edge 
of  the  sheet. 

While  it  is  believed  the  study  presented  here  provides  a basis  for 
tracing  certain  aspects  of  microcrack  growth  in  fibrous  composites, 
extensive  work  is  still  required  to  completely  understand  the  problem. 

A wider  range  of  fiber  flexibility  and  geometry  should  be  examined. 

The  nature  of  arrest  for  cracks  propagating  normal  to  the  fibers 
should  also  be  studied.  The  flexible  fiber  model  could  also  be  ex- 
tended by  a more  realistic  three-dimensional  model. 

To  predict  crack  growth  direction  in  the  three-dimensional  case, 
the  interaction  of  ail  three  crack  displacement  modes  will  have  to  be 
explored.  The  finite  element  method  may  prjve  to  be  promising  in 
this  problem;,  in  conjunction  with  the  total  energy  release  rate 
criterion.  This  criterion  itself  needs  further  experimental  verifi- 
cation, especially  with  materials  other  than  rigid  plastics. 
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APPENDIX 

Effective  Properties 


The  effective  properties  of  the  orthotropic  material  can  be  de- 
rived from  the  properties  of  the  constituent  materials  by  various  means. 
There  are  various  theories  proposed  that  are  based  on  descrete  element 
models,  variational  principles,  semi-empirical  approaches  and  others. 
These  methods  can  be  found  in  [3].  The  method  used  here  is  the  des- 
crete element  model.  Here,  each  constituent,  fiber  and  matrix,  is 
considered  separately  and  the  sum  of  their  individual  contributions  is 
assumed  to  be  the  same  as  that  of  the  orthotropic  material.  The  over- 
all stress-strain  relationship  used  here  is  for  an  orthotropic  body 
in  plane  stress  as  in  [142]. 


ei  = .£,  a.,  a.,  a..  = a..  , (i,j=l,2,6) 

1 3=1  13  3 13  31 


£,  = e,e=e,e,  = y 
1 x 2 y 6 xy 


(A— 1 ) 


cr,  = o , a = a , a,  = x 
1 x’  2 y 6 xy 


First,  an  element  of  the  material  is  considered  to  be  in  uniform 
strain  parallel  to  the  fibers.  One  considers  the  forces  in  each  fiber 
and  matrix  section,  and  sets  the  sum  of  all  the  forces  in  the  individual 
fiber  and  matrix  sections  equal  to  the  overall  force.  This  gives  the 
stiffness  in  the  fiber  direction  in  terms  of  the  volume  fractions  of 


the  constituents. 
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By  considering  the  lateral  contraction  of  the  overall  material,  one  can 
determine  that 


v10  = v V + v,  V. 
12  mm  f f 


(A-3) 


with  (A-l)  and  (A-2)  one  finds 


a12  V12  all 


(A-4) 


By  considering  an  element  of  the  material  in  simple  shear  parallel 
to  the  fibers,  the  shear  stresses  in  each  constituent  is  the  same. 

The  values  of  the  shear  displacement  in  each  fiber  and  matrix  section 
is  summed  and  set  equal  to  the  overall  shear  displacement.  In  terms 
of  the  fiber  and  matrix  shear  modulus  one  finds 


V,  V 
f , m 

a66  G.  + G 
f m 


(A-5) 


For  the  compliance  of  the  material  perpendicular  to  the  fibers, 
one  can  consider  an  element  of  the  material  to  be  loaded  by  an  overall 
force  perpendicular  to  the  fibers  and  with  no  overall  force  parallel 


to  the  fibers.  By  adding  the  displacements  perpendicular  to  the 
fibers,  for  each  fiber  and  matrix  section,  and  setting  the  sum  equal  to 


the  overall  displacement  of  the  element,  one  can  derive 


m . . 

a22  E 1 
m 


v (v  E - v E ) 
m m f f m 

E,  + EV 
f mm 


v.  (v  E - v E 
f f m m f 

E + E V 
m f f 


(A-6) 


These  compliances  were  used  as  the  effective 
continuum  composite  model  in  equation  (27). 
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